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We present MUSIC, an implementation of the Kurganov-Tadmor algorithm for relativistic 3+1 di- 
mensional fluid dynamics in heavy-ion collision scenarios. This Riemann-solver-free, second-order, 
high-resolution scheme is characterized by a very small numerical viscosity and its ability to treat 
shocks and discontinuities very well. We also incorporate a sophisticated algorithm for the de- 
termination of the freeze-out surface using a three dimensional triangulation of the hyper-surface. 
Implementing a recent lattice based equation of state, we compute pT-spectra and pseudorapidity 
distributions for Au+Au collisions at ^fs = 200 GeV and present results for the anisotropic flow 
coefficients V2 and 114 as a function of both pr and pseudorapidity 77. We were able to determine 
V4 with high numerical precision, finding that it does not strongly depend on the choice of initial 
condition or equation of state. 
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Hydrodynamics is perhaps the simplest description 
of the dynamics of a many-body system. Because it 
is coarse-grained, the complicated short-distance and 
short-time interactions of the particles are averaged out. 
Therefore, the effective degrees of freedom to describe 
the system reduce to a handful of conserved charges and 
their currents instead of some multiple of the number of 
particles in the system which can be prohibitively large. 
Yet, as long as the bulk behavior of a fluid is concerned, 
hydrodynamics is an indispensable and accurate tool. 

The equations of hydrodynamics arc thus simple: They 
are just the conservation laws and an additional equation 
of state (for dissipative hydrodynamics, constitutive re- 
lationships are also needed). In spite of their apparent 
simplicity, they can explain a vast amount of macroscopic 
physical phenomena ranging from the flow of water to 
the flight of airplanes. In this paper we are concerned 
in particular with applying ideal hydrodynamics to the 
description of extremely hot and extremely dense fluids 
- the Quark-Gluon Plasma (QGP) and hadron gasQ 

The idea that ideal hydrodynamics [l[ can describe the 
outcome of hadronic collisions has a long history start- 
ing from Landau Subsequent developments and 
applications to relativistic heavy-ion collisions have been 
carried out by many researchers [H-IHH and continue to 
this day. To describe the evolution of the system created 
by relativistic heavy-ion collisions, we need the following 
5 conservation equations 



ideal 



(e + V)u^u v - Vg^ v . 



« = 



(1) 
(2) 



where is the energy-momentum tensor and Jg is the 
net baryon current. In ideal hydrodynamics, these are 
usually re-expressed using the time-like flow 4-vector 



1 Wc will report on the extension of the current approach including 
viscous effects in another publication. 
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where e is the energy density, V is the pressure, ps is 
the baryon density and g^ v — diag(l, —1,-1,-1) is the 
metric tensor. The equations are then closed by adding 
the equilibrium equation of state 



V =V{s,p B ) 



(5) 



as a local constraint on the variables. 

In a first attempt to use these equations to study the 
QGP produced in relativistic heavy-ion collisions @ it 
was argued that at a very large y/s the boost invariant 
approximation should work well. Therefore one can elim- 
inate the longitudinal direction from the full 3D space. 
Further, it was assumed that the heavy ions are large 
enough so that the system is uniform in the transverse 
plane, thus eliminating all spatial dimensions from the 
equations. The energy-momentum conservation equation 
then simply becomes 



ds 



where t is defined as 



(6) 



together with the space-time rapidity rj s which trans- 
forms t, z coordinates to r, rj s coordinates as follows 



t = t cosh r/ s , 
z = t sinh r) s . 

If the equation of state is given by 

V = vie, 



(8) 



(9) 



where vl is the speed of sound squared, then we can easily 
find the solution of Eq.® 
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where Eq is the initial energy density at the initial time 
to- 

Although this solution is too simple to realistically de- 
scribe relativistic heavy-ion collisions, it still is a good 
first approximation for the mid-rapidity dynamics. How- 
ever, when one starts to ask more detailed questions 
about the dynamics of the evolving QGP such as the 
elliptic flow and HBT radii, it is not enough. One needs 
more sophisticated calculations. One of the first attempts 
to go beyond the Bjorken scenario was carried out in 
[56l H?]]. In the latter, the authors assumed cylindrical 
symmetry but otherwise used a fully three dimensional 
formulation and were successful in describing some SPS 
results available at the time. 

At SPS energy, the central plateau in the rapidity 
distribution is not very pronounced. It is more or less 
consistent with a Gaussian shape. In contrast, the cen- 
tral plateau extends over 4 units of rapidity at RHIC. 
Hence, as long as one is concerned only with the dynam- 
ics near the mid-rapidity region, boost invariance should 
be a valid approximation at RHIC, restricting the rele- 
vant spatial dimensions to the transverse plane. Pioneer- 
ing work on such 2+1 dimensional ideal hydrodynamics 
was carried out in pj, M, H, [H, S3 and gEHl- Much 
success has been achieved by these 2+1D calculations, in 
fact, too much to review in this work. Interested readers 
are referred to [34l |38| for a thorough review and exhaus- 
tive references. 

To go beyond 2+1D ideal hydrodynamics is chal- 
lenging. There are two main ingredients that need to 
be added - viscosities and the proper longitudinal dy- 
namics. Both require major changes in algorithm and 
computing resources. The main challenge for incor- 
porating viscosities into the algorithm is the appear- 
ance of the faster-than-light pr opagat ion of information. 
The Israel- Steward formalism 58|-6(3| avoids this super- 
luminal propagation [6lT[63j as does the more recent ap- 
proach in [30(. Since then a few groups have produced 
2+1D viscous hydrodynamic calculations. Two groups, 

0,IH|6l|6| and UlllillHi, use thc Isracl - Stewart for- 
malism of viscous hydrodynamics whereas another (48j 

uses the Ottinger-Grmela [66l - l68j formalism. There is 
much to discuss on the formalism of viscous hydrody- 
namics alone, but since it is not the main topic of this 
paper, we would like to defer the detailed discussion to 
our next publication where we will present our own vis- 
cous hydrodynamic calculations. 

The motivation to construct 3+1D hydrodynamics is 
to investigate the non-trivial longitudinal dynamics and 
its effects on the rapidity dependence of the transverse 
dynamics. Constructing a 3+1 dimensional ideal hydro- 
dynamics code, however, is not as simple as just adding 
one more dimension or one more equation to a code. 
The construction of a shock capturing algorithm, the 
freeze-out surface, all become much more intricate. So 
far there have been a few groups who have published 
their study of heavy-ion physics using realistic 3+1D 
ideal hydrodynamic simulations. One of them (29l. l32l |69| 



uses a fixed grid (Eulerian) algorithm to solve the hydro- 
dynamic equations, another |23l . l47j uses a Lagrangian 
approach which follows the evolution of each fluid cell. 
A somewhat different approach called smoothed particle 
hydrodynamics is used in Ref. [2(j ■ 

We have three major motivations to add another im- 
plementation to this list. The parameters such as the 
initial temperature profile and expansion rate can differ 
between the 2+1D calculations and the 3+1D calcula- 
tions and also among different approaches. Intuitively, it 
is clear that reality should favour 3+1D hydrodynamics. 
But since there are so many unknowns in the initial state, 
such as the exact initial energy density profile, initial flow 
profile and the initial baryon density profile, having an 
independent algorithm is important for verifying our un- 
derstanding of the initial condition and its uncertainty. 
This difference in 2+1D and 3+1D initial conditions is 
also important in jet quenching calculations since initial 
conditions can make a fair amount of difference in fixing 
the jet quenching parameters such as the initial temper- 
ature and the effective coupling constant. 

Another motivation is the desire to have a modular 
hydrodynamics code to which we can couple a high pr 
jet physics model such as in martini [70| . This is to 
examine the response of the medium to thc propagating 
jet as it loses energy to its surrounding medium. We are 
not yet at this stage but planning on implementing it in 
the near future. 

Last but not least, one should take advantage of recent 
progress in shock capturing algorithms to possibly sim- 
plify and certainly improve the calculations, creating an 
updated standard from which to assess the importance of 
viscous effects. The algorithm we use is usually referred 
to as Kurganov-Tadmor method [7l| . 

In this work we first review the Kurganov-Tadmor 
scheme (Section [H]) and present the implementation for 
relativistic ideal hydrodynamics in a three dimensional 
expanding geometry ( Section UTTj) . After discussing ini- 
tial conditions (Section IIVI) and the employed equations 
of state, which include a recent paramctrization of a com- 
bined lattice and hadron resonance gas equation of state 
(Section El, we introduce a new algorithm for determin- 
ing the freeze-out surface by discretizing the three dimen- 
sional hyper-surface into tetrahedra (Section lVij) . Finally 
we show first results for particle spectra including reso- 
nance decays from resonances up to 2 GeV, elliptic flow, 
and the anisotropic flow coefficient V4 (Section IVII[) . It 
is demonstrated that the latter is highly sensitive to dis- 
cretization errors which are shown to be well under con- 
trol for fine enough lattices. 



II. KURGANOV-TADMOR METHOD 

Hydrodynamic equations stem from conservation laws. 
Hence, they take the following general form: 



dtp a 



-V-J 0! 



(11) 



3 



where a runs from to 4, labelling the energy, 3 com- 
ponents of the momentum and the net baryon density. 
The task is then to solve these equations together with 
the equation of state. Even though they have a decep- 
tively simple form, they are remarkably subtle to solve. 
In this section, we briefly sketch the Kurganov and Tad- 
mor scheme (KT) [7l[, which we use for the solution of 
Eq.(HU). 

To illustrate the method, consider the following single 
component conservation equation in 1 spatial dimension 



d t p = -d x J, 



112) 



together with an equation that relates J to p such as 
J = up. All the essential features of KT can be ex- 
plained with this simple example. As it was shown in 
Ref. [7l| , higher dimensions can be dealt with by simply 
repeating the treatment here for all spatial dimensions. 
Coupled conservation equations make the calculation of 
the maximum local propagation speed more complicated 
(see below), but there is no conceptual complication in 
doing so. 

The need for more sophisticated numerical methods 
in solving conservative equations in part comes from the 
fact that a naive discretization of Eq. (|T2"j) such as 



pT 



At 



jn jn 

J j+i J j-i 

2Ax 



(13) 



with J = vp is unconditionally unstable. That is, the 
solution will cither grow without bound as t increases or 
start to oscillate uncontrollably. Here the superscript n 
indicates that the quantity is evaluated at t n = to + nAt 
and the subscript j indicates that the quantity represents 



the value at 



devises a scheme where numerical damping is introduced. 
For instance, suppose one replaces p™ in the left hand side 
of Eq. (fT3]) with the spatial average (p™ + i + p"_i)/2. In 
the small At and Ax limit, this well-known Lax method 
[72l |73| is equivalent to solving 



d t p = -d x J 



[Axy 
2At 



d x p- 



(14) 



The second term is the numerical dissipation term of- 
ten referred to as the "numerical viscosity". Different 
schemes introduce different forms of the numerical vis- 
cosity term. 

This simple method does stabilize the numerical solu- 
tions but one also can immediately see that (Ax) 2 /At 
must not be large. Otherwise, this artificial term will 
dominate the numerical evolution of the system. There- 
fore in this method, a finer time resolution result can- 
not be computed without making the number of spatial 
grid points correspondingly large. Many other numerical 
methods also have a 1/ At behavior for the artificial vis- 
cosity, including KT's immediate predecessor J74|. How- 
ever, in KT the artificial viscosity does not depend on 
At. It only depends on some positive power of Ax and 
we are free to take the At — >• limit. As a bonus, this 



allows us to turn a set of difference equations into a set of 
ordinary differential equations as explained below. This 
places the vast array of ODE solvers at one's disposal, 
thus making this method much more versatile. 

Notably, KT is a MUSCL-type (Monotonic Upstream- 
centered Schemes for Conservation Laws) finite volume 
method (75j in which the cell average of the density p 
around Xj is used instead of the value of the density at 
Xj . Then, the conservation equation for the cell average 



X 3 + l/2 



dx p(x, t) , 



(15) 



becomes 



d n M _ J(xj-i/2,t) - J(x j+1/2 ,t) 
dt PA '~ Ax ' 1 ' 

and the current and charge density at values other than 
the Xj are contructed using a pieccwisc linear approxima- 
tion. This method leads to discontinuities at the halfway 
points Xj±i/2 where the current is evaluated. Kurganov 
and Tadmor solved this problem using the maximal lo- 
cal propagation speed a — \dJ/dp\ to identify how far 
the influence of the discontinuities at Xj±i/2 could travel, 
and divided the space into two groups; one with elements 
that include a discontinuity and one where the solution is 
smooth. The exact procedure of doing this is explained 
in Appendix [A"l 

Here, we quote Kurganov and Tadmor's final result for 
the conservation equation in the At — >• limit: 



dt 



Pj(t) 



#7 + 1/2 ft) - .Hj-l/gft) 

Ax 



(17) 



jAx. One can make this stable if one where 



H 



i±i/2 



J {Xj±l/2,+ ,t) + J (Xj±l/2- 1 t) 



Oj±l/2(i) 



(pj±l/2,+ (t)-Pj±l/2,-(*)) ,(18) 



with 



Pj+i/2,+ - Pj+i 2~(Px)j+i 

Ax 

Pi+l/2- = Pi + —{Px)j ■ 



(19) 
(20) 



The order of the spatial derivatives {p x )j is chosen by 
the minmod flux limiter 



(p x )j = minmod 



where 



J'.- • 1 _ Pi Pj+i - Pi-i a Pi ~ ft-i 



Ax 



2Ax 



Ax 



minj {xj}, if Xj > Vj 
minmod(xi, X2, • • • ) = ^ maxjjxj}, if Xj < Vj 

0, otherwise 

and 1 < 8 < 2 is a parameter that controls the amount 
of diffusion and the oscillatory behavior. This is also 
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our choice with 8 = 1.1. This allows for higher accu- 
racy using the second-order approximation where possi- 
ble and avoids spurious oscillations around stiff gradients 
by switching to the first order approximation where nec- 
essary. 

The Kurganov-Tadmor method, combined with a suit- 
able flux limitcr such as the one just described, is 
a non-oscillatory and simple central difference scheme 
with a small artificial viscosity which can also handle 
shocks very well (for an extensive comparison with other 
schemes in this regard, see Ref. [Z6|). It is also Riemann- 
solver free and hence does not require calculating the 
local characteristics. This scheme is ideally suited for 
hydrodynamics studies. 



III. IMPLEMENTATION 

We now describe our implementation of the KT algo- 
rithm for rclativistic heavy-ion collisions, dubbed MUSIC, 
MUScl for Ion Collisions. 

As in most ideal hydrodynamics implementations for 
heavy-ion collisions, the most natural coordinate system 
for us is the r — rj s coordinate system defined by Eq. (J5]) . 
In the T—r] s coordinate system, the conservation equation 

J M = becomes 



where 



d T {Tj T )+d Vs J^ +d v {Tj V ) =0, 

J T = (cosh J° — sinh?7 s J 3 ) , 
J' h = (cosh J 3 - sinh77s J°) , 



(21) 



(22) 
(23) 



which is nothing but a Lorentz boost with the space- 
time rapidity rj s = tanh _1 (z//j). The index v and w in 
this section always refer to the transverse x, y coordinates 
which are not affected by the boost. Applying the same 
transformation to both indices of T Ml/ , one obtains 



d T (TT TT ) + Vs (T"= T ) + O v (tT vt ) + T"^ = , (24) 



and 



d T {TT T 'i°) + Vs (T"'"') + d v {rT v ^) + T T '" S = ,(25) 



d T (TT TV ) + d rh (T nsV ) + d w (rT wv ) = , (26) 



and 



These 5 equations, namely Eq. (|2Tj) for the net baryon 
current, and Eqs. (|2~4l |2~51 |2"BT) for the energy and momen- 
tum are the equations we solve with the KT scheme ex- 
plained in the previous section. Multi-dimension is dealt 
with by repeating the KT scheme in each direction [7l| . 
The source term is dealt with by follo wing the suggestions 
in the original KT paper and others 771 ] . 

At each time step, the new values of J T ,T TT , T TTIs and 
T TV are obtained by solving the semi-discrete version 
of KT using Heun's rule. Heun's rule is a form of the 



second-order Rungc-Kutta method which can be stated 
as follows. Suppose we have a differential equation 



(27) 



A numerical solution of this equation can be obtained by 
applying the following rules 

1. Compute ki = f(t,p n ). 

2. Compute p' n+ i = p n + &iAt. 

3. Compute fc 2 = f(t + At,p' n+1 ). 

4. Compute p n+1 = p n + (hi + k 2 )At/2. 

Once new values of J r , T TT , T ,,sT and T VT are obtained, 
the following ideal gas expressions 



T TT = {e + V)u T u T - V 
T Trh = {e + V)v?*u T 
T TV = (e + V)u T u v 

,r = P u T 

together with the equation of state 
V = V(e,p) 



(28) 
(29) 
(30) 
(31) 

(32) 



determine the net baryon density p, the pressure V, the 
energy density e, and the flow velocity u^. The flow com- 
ponents u T and u r?s here are given by the Lorentz boost 
with the space time rapidity ?7s exactly as in Eqs. (|221 
|2"3"|) . Hence, they still satisfy the normalization condition 

U 2 T = l + U 2 x +U 2 y + (l/T 2 )u 2 nB . 

Explicitly, the values of e and p are obtained by itera- 
tivcly solving the following coupled equations 



e = T TT - 



P=J T \ 



1 e + V{e lP ) 
T TT + V(e,p) ' 



(33) 
(34) 



where K = (T^ T ) 2 + (T^f + (T« r ) 2 . A good initial 
guess turned out to be either the value at the previous 
time step or just the initial T TT and J T . Knowing e, p, we 
can calculate the pressure V = V(e,p) and u T = J T /p. 
These then determine the spatial flow vector components 
as 



(35) 



(e + T)u T 



for i = rj s , x, y. With these e, p, V and m m , the whole 
can be reconstructed and be used at the next time step. 

In addition to the currents, we need to find the maxi- 
mum local propagation speed at each time step. The 
maximum speed in the k direction is given by the maxi- 
mum eigenvalue of the following Jacobian: 



a jk 
qk _ UJ a 
Jab gjr 



(36) 
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where J£ with a = 0, 1, ■ • • ,4 stand for the 5 currents 
(net baryon, energy and momentum) . The whole matrix 
is quite complicated. However, with the help of math- 
ematica [731, ^ turned out that the eigenvalues can be 
analytically calculated. 

If there is no net baryon current to consider, two of 
the 4 eigenvalues in the k = x,y direction are degenerate 
and equal to u k ju T . The remaining two are 



A; 



a±Vb 



D 



(37) 



with 



.4 
B 
D 



u T u k (l-v 2 s ), 



— u h — I U T u. 



(u 2 T - 1) v\ 



1) <} v 



2 

S 5 



(38) 

where v 2 = T"{e) is the speed of sound squared. In the 
k = rj s direction, we have the same expression, but the 



eigenvalues are scaled with 1/r, that is A,, 



A 



It. 



The same expressions for the Cartesian case was obtained 
in [79|. The largest eigenvalue is thus 



IA 



largest I 



D 



D 



(39) 



with an additional 1 /r factor for k = i] s . 

Even with the net baryon current present, the eigen- 
values of the the resulting 5x5 matrix can be computed 
analytically. Consider first k = x,y directions. Among 
the 5 eigenvalues, 3 are degenerate and equal to u k /u T . 
The other two are again given by Eq. (j3"8"|) with 

v 2 = d e T + (p/(e + T))d p T. (40) 

It is obvious that the p — >• limit coincides with the 
no current case. In the t] s direction, the expressions for 
the eigenvalues are the same except for the overall scale 
factor 1/r. The maximum eigenvalue is again given by 
Eq. (f3l)| with the above substitution. 

Importantly, MUSIC is fully parallelized to run on many 
processors simultaneously. To achieve this, the lattice is 
truncated in the rj s direction so that each processor only 
has to evolve the system on one slice, communicating the 
cell values at the boundary to the neighboring processors 
every time step. This leads to an increase in speed almost 
(minus the time necessary for communication between 
the processors) linear in the number of processors used. 

The typical size of a time step is of the order of 
0.01 fm/c. Energy conservation is fulfilled to better than 
1 part in 30,000 per time step. 



IV. INITIAL CONDITIONS 

The initialization of the energy density is done using 
the Glauber model (see [8(| and references therein): Be- 
fore the collision the density distribution of the two nuclei 
is described by a Woods-Saxon parametrization 

Po 



l + exp[(r- R)/d] 



(41) 



with R = 6.38 fm and d — 0.535 fm for Au nuclei. The 
normalization factor po is set to fulfill J d 3 rp A (r) = A. 
With the above parameters we get po = 0.17 fm~ 3 . The 
relevant quantity for the following considerations is the 
nuclear thickness function 



T A {x,y) 



dz p A (x,y,z) 



(42) 



where r = \J x 2 + y 2 + z 2 . The opacity of the nucleus 
is obtained by multiplying the thickness function with 
the total inelastic cross-section er of a nucleus-nucleus 
collision. 

Experiments at SPS found that the number of final 
state particles scales with the number of wounded nucle- 
ons, nucleons that interact at least once in the collision. 
Deviations from the scaling are observed at RHIC. 

Statistical considerations allow to express the number 
of wounded nucleons in the transverse plane by the nu- 
clear thickness function of one nucleus, multiplied with a 
combinatorial factor involving the nuclear thickness func- 
tion of its collision partner. This factor ensures that the 
participating nucleon does not penetrate the finite oppos- 
ing nuclear matter without interaction. For noncentral 
collisions of nuclei with mass numbers A and B at im- 
pact parameter 6, the number of wounded nucleons per 
transverse area is given by (27| 



nwN(x,y, 6) 



T A (x+-,y) 



1 _f 1 _ a n T B{x -\,y) 



T B (x- -,y) 



a T A (x + |,y) 



(43) 



Integrating this expression over the transverse plane 
yields the total number of wounded nucleons (partici- 
pants) as a function of the impact parameter. We com- 
pute the relevant quantities using routines adapted from 

LEXUS [IH. 

At high energies the density of binary collisions be- 
comes of interest. After suffering their first collision, the 
partons travel on through the nuclear medium and are 
eligible for further (hard) collisions with other partons. 
This leads to the notion that one has to count the binary 
collisions. The density of their occurence in the trans- 
verse plane is simply expressed by the product of the 
thickness function of one nucleus with the encountered 
opacity of the other nucleus, leading to 

n BC {x, y, b) = <r T A {x + 6/2, y)T B {x - 6/2, y) . (44) 

The total number of binary collisions shows a stronger de- 
pendence on the impact parameter than does the number 
of wounded nucleons. 

We now assume that the initial state of matter in the 
transverse plane is governed entirely by the physics of 
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'soft' and 'hard' processes represented in terms of the 
densities of wounded nucleons and binary collisions, re- 
spectively. Shadowing effects by the spectators do not 
play a role at RHIC energies because the spectators leave 
the transverse plane at z = on a timescale of less than 
1 fm/c. 

Whether the deposited energy density or entropy den- 
sity scales with the density of wounded nucleons or binary 
collisions is not clear from first principles. As mentioned 
above, SPS data suggests that the final state particle 
multiplicity is proportional to the number of wounded 
nucleons. At RHIC energies a violation of this scaling 
was found (the particle production per wounded nucleon 
is a function increasing with centrality. This is attributed 
to a significant contribution from hard processes, scaling 
with the number of binary collisions). 

We parametrize the shape of the initial energy density 
distribution in the transverse plane as 

W(x, y, b) = (1 - a) r7, WN (a;, y, b) + a n BC {x, y, b) , (45) 

where a determines the fraction of the contribution from 
binary collisions. 

Alternatively, we can scale the entropy density as op- 
posed to the energy density as in (|4"5")) . This leads to more 
pronounced maxima in the energy density distributions 
because of the relation e ~ s 4 / 3 in the QGP phase. How- 
ever, a similar effect can be achieved by increasing the 
contribution of binary collision scaling a. 

For the longitudinal profile we employ the prescription 
used in m iHl El ii, US]- It is composed of two 
parts, a flat region around r/ s — and half a Gaussian in 
the forward and backward direction: 



exp 



2o* 



M - JJflat/2) 



(46) 



The full energy density distribution is then given by 

e(x, y, Vs ,b)= £ H( Vs ) W(x, y, b)/W(0, 0, 0) . (47) 

The parameters r/fi at and a v are tuned to data and will 
be quoted below. 

V. EQUATION OF STATE 

To close the set of equations (jUJ [H [231 I25J) wc must 
provide a nuclear equation of state V(e,p) which relates 
the local thermodynamic quantities. We present calcu- 
lations using a modeled equation of state (EOS-Q) also 
used in AZHYDRO [HI, [33, IH] as well as one extracted 
from recent lattice QCD calculations [Hj]. 

For the EOS-Q, the low temperature regime is de- 
scribed as a non-interacting gas of hadronic resonances, 
summing over all resonance states up to 2GeV f87| . 
Above the critical temperature T cr it = 164 MeV, the sys- 
tem is modeled as a non-interacting gas of massless u, d, s 
quarks and gluons, subject to an external bag pressure B. 



The two regimes are matched by a Maxwell construction, 
adjusting the bag constant i? 1 / 4 = 230 MeV such that for 
a system with zero net baryon density the transition tem- 
perature coincides with lattice QCD results [III [89| . The 
Maxwell construction inevitably leads to a strong first 
order transition, with a large latent heat. 

However, lattice results suggest a smoother transition. 
Recently, in [86j several parametrizations of the equa- 
tion of state which interpolate between the lattice data 
at high temperature and a hadron-resonance gas in the 
low temperature region were constructed. We adopt the 
parametrization "s95p-vl" (and call it EOS-L in the fol- 
lowing), where the fit to the lattice data was done above 
T = 250 MeV, and the entropy density was constrained 
at T = 800 MeV to be 95% of the Stefan-Boltzmann 
value. Furthermore, one "datapoint" was added to the fit 
to make the peak in the trace anomaly higher. See [86j 
for more details on this parametrization of the nuclear 
equation of state. 



VI. FREEZE-OUT 

The spectrum of produced hadrons of species i with 
degeneracy gi is given by the Cooper- Frye formula (90| : 

E 7jT = a d f 72 = 9i [ , (48) 

d A p dyp T dp T d(pp J s 

with the distribution function 

1 1 



/(«%) 



(2tt) 3 exp((u"p„ - hJ/Tfo) ± 1 



(49) 



We assumed that at freeze-out every infinitesimal part 
of the hyper-surface E behaves like a simple black body 
source of particles (this assumption will be modified when 
including viscosity). The collective velocity of the fluid 
on the hyper-surface, which results from longitudinal and 
transverse flow, is taken into account by using the invari- 
ant expression E = E(x) = u fJ, (x)p f ^. To evaluate the 
right hand side of (|38]l we need to determine the freeze- 
out hyper-surface 

E = (Y, (x,y,ri s ),j: 1 (x,y,r] s ),j: 2 (x,y,r] s ),j: 3 (x,y,'n s )) 
= ( T f( x , V, Vs)coshT] s , x, y, T f (x, y, ii s )sm\nj s ) , (50) 

where Tf{x,y,rj s ) is the freeze-out time, determined by 
when the energy density (or temperature) falls below the 
critical value £fo (or Tpo)- The normal vector on this 
surface is given by 



c^E, 



dY, v <9E dYP 



dxdydrjs , (51) 



dx dy drj s 

with the totally anti-symmetric tensor of fourth order 



1 even permutation 
i^uXp = <. —1 odd permutation (52) 
otherwise 
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Using (|50j) we find 

S " = U: Smh% + TS co ^,-r f — , -r f — , 



dr 



— — cosh n s — Tr sinh ri s dxdydn 
oris 



(53) 



To evaluate Eq. (|48|) we need to determine and 
p^dSp. The hydrodynamic evolution calculation pro- 



vides u T ,u x ,u v , u 71 * , so we express p^ ai 



Tot cosh(y — ?7 S 
sinn(y — rj s 



with tot = \J m 2 + pi?, where to is the mass of the con- 
sidered particle, and obtain 



u^Pfj, = u T p T 



u x p x 



u y p y _ T i u is p vs 



u 1 tut cosh(y — r/. 
— tu Vs uit sinh(y — r) s ) 



u x p x 



u v p y 



(54) 



Note that y's appearing in the cosh and sinh functions 
represent the rapidity of the produced hadron. We can 
express £ from Eq. (|50|) in terms of r — rj s coordinates 

S Q = 5^, 5>) = (T f (x,y,r, s ),x,y,r, 3 ), (55) 

and get 

d 3 S Q = (1,— -i, — , ) V det q dxdydrjs 
ox ay orj s 

with g the metric in r — rj s coordinates. The scalar prod- 
uct of (f5"6) with p Q is then found to be 



P a d^ c 



m Tj— fa sinh(y - rj s )) - r f p T ■ V t t/ 



dxdydr] s , 
(57) 



with the two-dimensional derivative Vt = (d x ,d y ). 

In the limit that 77 does not depend on rj s we recover 
the Bjorken result 

p a d i Y* a = ^totcos1i(j/ — r) s ) — pr ■ VTTfjTfdxdydr] s . 

(58) 

In practice, we need to determine d 3 E Q geometrically. 
In previous works a simple algorithm has been used 
}32l HlJ that adds a cuboidal volume element to the to- 
tal freeze-out surface whenever the surface crosses a cell, 



e.g., if the quantity e — epo changes sign when moving 
along the x direction, one adds a volume element of size 
Aj/A?/ s At with its hyper-surface vector pointing in the 
x direction (towards lower energy density). So surface 
vectors always point along one of the axes x, y, rj s , or 
t. This method overestimates the freeze-out surface it- 
self but is sufficient for computing particle spectra. How- 
ever, it turns out that for computing anisotropic flow and 
especially higher harmonics than V2 it is essential to de- 
termine the freeze-out surface much more precisely. To 
do so, within MUSIC we employ the following method: 

We define a cube in 4 dimensions that may reach over 
several lattice cells in every direction and over several 
r-stcps, and determine if and on which of the cube's 32 
edges the freeze-out surface crosses. In this work we let 
the cube extend over one lattice cell in each spatial di- 
mension and over 10 steps in the time direction. If the 
freeze-out surface crosses this cube, we use the intersec- 
tion points to perform a 3D-triangulation of the three di- 
mensional surface element embedded in four dimensional 
space. This leads to a group of tetrahedra, each con- 
tributing a part to the hyper-surface-vector. This part is 
of the form 



d^ = e 



A <x B P C -f/Q 



(59) 



where A, B, and C are the three vectors that span the 
tetrahedron n. The factor 1/6 normalizes the length of 
the vector to the volume of the tetrahedron. We demand 
that the resulting vector points into the direction of lower 
energy density, i.e., outwards. The vector-sum of the 
found tetrahedra determines the full surface- vector in the 
given hyper-cube. 

Depending on where the freeze-out surface crosses the 
edges, the structure may be fairly simple (e.g. 8 crosses, 
all on edges in x-direction) or rather involved (crossings 
on edges in many different directions). The current al- 
gorithm is close to perfect and fails to construct hyper- 
surface elements only in very rare cases. Typically these 
are cases when the surface crosses the cube in many dif- 
ferent directions, e.g. in the rj s , x, and r direction. How- 
ever, even for these cases a full reconstruction can usually 
be achieved and the algorithm was found to succeed in 
determining the volume element in ~ 99% of the cases 
for the studied systems. The ~ 1% of surface elements 
that could not be fully reconstructed usually miss only 
one tetrahedron. Because one typocally needs between 8 
and 20 tetrahedra to reconstruct a cell, the error intro- 
duced by missing one tetrahedron in the 1% of the cells 
lies between 5 and 15%. Considering the high complexity 
of the triangulation procedure in four dimensions, this is 
a very satisfactory result. 



VII. RESULTS 



2 In this section, our definition of the 4-vector component v Vs 
(coshes v s — sinh»7 s v°)/t carries an extra factor of 1/r. 



To obtain results for particle spectra, we first compute 
the thermal spectra of all particles and resonances up to 
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FIG. 1. (Color online) pr spectra for n , K , and p at central 
collisions using different equations of state (thin lines: AuAu- 
I (EOS-Q), thick lines: AuAu-3 (EOS-L)) compared to 0-5% 
central PHENIX data [95|. The used impact parameter was 
b = 2 Aim. 



FIG. 2. (Color online) Centrality dependence of pseudorapid- 
ity distribution compared to PHOBOS data (97|. From top to 
bottom, the used average impact parameters are b — 2.4 fm, 
b = 4.83 fm, b = 6.7, fm, and b = 8.22 fm. 



~ 2 GeV using Eq. (j4*8"f and then perform resonance de- 
cays using routines from AZHYDRO [2l|, [H, |H, |93| that we 
generalized to three dimensions. Unless indicated other- 
wise, all shown results include the resonance feed-down. 
Typically, the used time step size is At sa 0.01 fm/c, and 
the spatial grid spacings are Aa; = Ay = 0.08 fm, and 
Ar/ S = 0.3. This is significantly finer than in previous 
94 1 for example uses Ar = 0.3 fm/c, 
0.3. The possibility 



3+1D simulations 
Ax = Ay = 0.3 fm, and A77. 



to use such fine lattices is an improvement because it is 
mandatory when computing higher harmonics like V4 as 
demonstrated below. Another advantage of using large 
lattices is that in the KT scheme the numerical viscos- 
ity decreases with increasingly fine lattices (sec Appendix 
[A"]) . The spatial extend of the lattice used in the follow- 
ing calcualtions is 20 fm in the x and y direction, and 20 
units of rapidity in the r) s direction. 



A. Particle spectra 

In Fig.[T] we present the transverse momentum spec- 
tra for identified particles in Au+Au collisions at y/s = 
200 GeV compared to data from PHENIX [H]. The used 
parameters are indicated in Table [J] They were obtained 
by fitting the data at most central collisions. 

We reproduce both pion and kaon spectra well. The 
model assumption of chemical equilibrium to very low 
temperatures leads to an underestimation of the anti- 
proton spectrum. The overall shape is however well re- 
produced, even more so with the EOS-L that leads to 
flatter spectra [86[. 

One way to improve the normalization of the proton 
and anti-proton spectra (as well as those of multistrangc 
baryons) is to employ the partial chemical equilibrium 



model (PCE) [32|, Hfi, [96|, which introduces a chemi- 
cal potential below a hadron species dependent chemical 
freeze-out temperature. Note that the initial time was 
set to To = 0.4 fm/c when using the EOS-L to match the 
data. The quoted parameter sets fit the data very well, 
however, they do not necessarily represent the only way 
to reproduce the data and a more detailed anaylsis of 
the whole parameter space may find other parameters to 
work just as well. 

Next, we show the pseudorapidity distribution of 
charged particles at different centralities compared to 
PHOBOS data [13 in Fig.d The only parameter that 
changes in going to larger centrality classes is the im- 
pact parameter. Experimental data is well reproduced 
also for semi-central collisions, showing that the results 
mostly depend on the collision geometry. The used im- 
pact parameters, b = 2.4 fm, b = 4.83 fm, b = 6.7, fm, 
and b = 8.22 fm, were obtained using the optical Glauber 
model and correspond to the centrality classes used by 
PHOBOS. We show the centrality dependence of the 
transverse momentum spectrum of 7r~ in Fig. [3] Devi- 
ations occur for more peripheral collisions because the 
soft collective physics described by hydrodynamics be- 
comes less important compared to jet physics in peri- 
pheral events. However, we find smaller deviations than 
e.g. H3. 



Finally we present results for the average transverse 
momentum of pions and kaons as a function of pseudo- 
rapidity in central collisions. We compare with — 5% 
central data by BRAHMS [98| and find good agreement 
for kaons, but slightly larger values for pions. This could 
be expected because the calculated pr spectra are slightly 
harder than the experimental data, especially when using 
the EOS-L (sec Fig.H]). 
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EoS 


to [fm] 


eo[GeV/fm 3 ] 


po[l/fm 3 ] 


e FO [GeV/fm 3 ] 


T F o[MeV] 


Of 




a v 


AuAu-1 


EOS-Q 


0.55 


41 


0.15 


0.09 


« 130 


0.25 


5.9 


0.4 


AuAu-2 


EOS-Q 


0.55 


35 


0.15 


0.09 


« 130 


0.05 


6.0 


0.3 


AuAu-3 


EOS-L 


0.4 


55 


0.15 


0.12 


« 137 


0.05 


5.9 


0.4 



TABLE I. Parameter sets. 




FIG. 3. (Color online) Centrality dependence of ir~ transverse 
momentum spectra compared to PHENIX data [95| . The 
curves (both data and hydro) for 10 — 15%, 15 — 20% and 
20 — 30 % centrality are scaled by a factor of 5, 25, and 150, 
respectively. Thick lines are for parameter set AuAu-3 (EOS- 
L), thin lines for AuAu-1 (EOS-Q). 



FIG. 4. (Color online) (pt) for positive kaons and pions as a 
funciton of rapidity compared to most central BRAHMS data 
. The used impact parameter is 6 = 2.4 fm. Different lines 
correspond to different parameter sets: From top to bottom: 
AuAu-3 (EOS-L), AuAu-1, AuAu-2 (EOS-Q). 



B. Elliptic flow 

Wc present results for i; 2 as a function of px integrated 
over the pseudorapidity range —1.3 < r\ < 1.3, which cor- 
responds to the cut in the analysis by STAR [99| that we 
compare to. We show results for identified hadrons ob- 
tained using parameter set AuAu-1 (EOS-Q) and AuAu- 
3 (EOS-L) in Fig.|5] While the pion elliptic flow is rela- 
tively well described for both equations of state, we find 
an overestimation of the anti-proton v-i , especially when 
using the EOS-L. This is compatible with results in [86| . 

Charged hadron V2 is presented in Fig.|6]where we com- 
pare results using different contributions of binary colli- 
sion scaling a which lead to different initial eccentricities. 
We also show the result obtained by using the EOS-L, 
which is somewhat above the EOS-Q result for lower px 
but bends more strongly to be smaller at pr = 2 GeV. 

Overall, we find that while the pion V2 is well repro- 
duced, both anti-proton and charged hadron w 2 is over- 
estimated for both parameter sets. So there is room for 
viscous corrections that have been found to reduce at 
p T = 1.5 GeV by 20% for ?7 shoar /s = 0.08 jH H3, HI Hi ■ 

Fig. [7] shows V2 of positive pions for different central- 
ity classes, again comparing calculations using parameter 
sets AuAu-1 and AuAu-3 with experimental data [99| . 
In both cases the agreement with the experimental data 
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FIG. 5. (Color online) pt dependence of the elliptic flow 
coefficient V2 for tt~ and p using parameter set AuAu-1 (EOS- 
Q, thin lines) and AuAu-3 (EOS-L, thick lines) compared to 
STAR data from [g|. 



that is available to up to pt = 1 GeV is very reasonable. 

In Fig. [8] we present the result for V2 as a function 
of ps eudorapidity ry, comparing to data from PHOBOS 
jlOOj . As earlier calculations [H, |47[ the hydrodynamic 
model calculation overestimates the elliptic flow espe- 
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FIG. 6. (Color online) pr dependence of the elliptic flow co- 
efficient V2 for charged hadrons using parameter sets AuAu-1, 
AuAu-2 (EOS-Q), and AuAu-3 (EOS-L) compared to STAR 
data from [99l ]. 



FIG. 8. (Color online) Pseudorapidity dependence of the el- 
liptic flow coefficient Vi for charged hadrons using parameter 
sets AuAu-1 (EOS-Q), and AuAu-3 (EOS-L) compared to 
PHOBOS data from [100|. 



C. Higher harmonics 
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FIG. 7. (Color online) Elliptic flow coefficient V2 for positive 
pions using parameter set AuAu-1 (EOS-Q, thin lines) and 
AuAu-3 (EOS-L, thick lines) for different centrality classes 
compared to STAR data from [gjj. 



cially at forward and backward rapidities. This is most 
likely due to the fact that the assumption of ideal fluid 
behavior is no longer valid far away from the midrapidity 
region. Calculations combining hydrodynamic evolution 
with a hadronic after-burner improve on this [111 l47j . 
Effects of viscosity on ^2(77) in the 3+1 dimensional sim- 
ulation have been estimated to be stronger at larger \r/\ 
[l0ll | and it will be interesting to see what a full compu- 
tation will yield. 



The extraction of higher harmonic coefficients from the 
computed particle distributions has to be done with great 
care. Apart f rom being highly sensitive to the initial 
conditions Il02| , the fourth harmonic coefficient V4 is 
also highly sensitive to the discretization of the freeze- 
out surface and lattice artifacts. Where other quantities 
such as pt spectra and V2 are almost unaffected by a 
change of the lattice resolution or the freeze-out method, 
V4 depends strongly on the method and the lattice spac- 
ing. Using the simplified freeze-out surface algorithm de- 
scribed above, the dependence of V4 on the discretization 
becomes very strong (in this case V4 is negative when us- 
ing a 128 3 -lattice and only becomes positive and slowly 
approaches the correct value for much finer lattices). 

It is therefore necessary to work on very fine lattices 
and have a very sophisticated algorithm for determin- 
ing the freeze-out surface in order to obtain reliable re- 
sults for 1*4. To measure the error introduced by the 
anisotropic discretization of the lattice (lattice along the 
diagonal in the transverse plane looks different than along 
one of the axes), we compute twice: once with the im- 
pact paramter along the x-axis, once with the impact 
parameter along the diagonal in the rr-y-planc. The dif- 
ference between the results is a measure of discretization 
errors in V4 and is shown for the pion V2 in Fig.[SJ The 
difference decreases significantly when going from a 64 2 
to a 320 2 lattice in the transverse plane. Hence, the nu- 
merical error of U4 is well under control. 

Fig-Hni shows V4 of charged hadrons computed with 
both parameter set AuAu-1 (EOS-Q) and AuAu-3 (EOS- 
L). We added error bands representing an estimate for 
the discretization error on the used 256 2 x 64 lattice. Mo- 
tivated by the results shown in Fig. HO we choose ±15%. 
Experimental data for mid-central centrality classes is 
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FIG. 9. (Color online) Pion-i>4 (without resonance decays) as 
a function of pr computed on different lattices. The upper 
curve of each band is the result from when the impact param- 
eter is aligned with the diagonal in the x-y-pl&ne, the lower 
curve from when it is aligned with the a:-axis. The absolute 
value of the impact parameter is b = 6.3 fm. We compare to 
charged hadron STAR data from [9!| ■ 



FIG. 11. (Color online) V4 for charged hadrons using parame- 
ter sets Au Au-1 and AuAu-3 at 6 = 7.5 fm compared to STAR 
data from Il03ll. 



mostly from our overestimation of «2 at high pr, while 
t>4 is well reproduced (see Figs. 151 and rTU|) . 
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FIG. 10. (Color online) Vi for charged hadrons using all pa- 
rameter sets for different centrality classes compared to STAR 
data from 19911. See text for details. 



well r epro duced in both cases, and contrary to expecta- 
tions jl02| , we find that 1*4 is not very sensitive to either 
the initial condition or the equation of state. This is 
also visible in Fig.[TT] where we show V4 as a function 
of ps eudorapidity compared to preliminary STAR data 
[1031 ] . Also here we add an error band to indicate the 
discretization error, estimating it to be ±15%. 

We have checked that the ratio V4/v% approaches 0.5 
for large pr as it should for an id eal fl uid, at least in 
the limit of small impact parameter jl04| . The difference 
to the data, which for charged hadrons in mi nimu m bias 
collisions is about constant at 1*4/1*2 « 1.2 99lll05l ] comes 



VIII. CONCLUSIONS AND OUTLOOK 

We presented first results from our newly developed 
3+1 dimensional rclativistic fluid dynamic simulation, 
MUSIC, using the Kurganov-Tadmor high-resolution cen- 
tral scheme to solve the hydrodynamic equations. The 
method handles large gradients very well, which makes it 
ideal for future explorations of 'lumpy' initial conditions 
or energy-momentum deposition by jets. It also has a 
very small numerical viscosity, which is a prerequisite for 
extracting physical viscosities in the future extension to 
dissipative hydrodynamics. 

We showed a detailed comparison of results using 
different equations of state including a very recent 
paramctrization of combined lattice and hardon reso- 
nance gas equations of state. Our calculations of identi- 
fied hadron py-spectra, pseudorapidity distributions and 
elliptic flow coefficients of charged hadrons in Au+Au 
collisions at the highest RHIC energies reproduced results 
of earlier 3+1 dimensional simulations. In addition, we 
were able to obtain reliable results for the anisotropic flow 
coefficient v^, which is highly sensitive to discretization 
errors. This was possible by developing a sophisticated 
algorithm for determining the freeze-out surface in four 
dimensions and running the simulation on fine lattices 
on many processors parallely to obtain results within a 
reasonable amount of time. We found that contrary to 
earlier expectations, is not very sensitive to the initial 
conditions. Neither is it very sensitive to the equation of 
state. 

The next step will be the inclusion of viscous effects 
which we will present in a forthcoming work. We also 
plan to combine the simulation with our event generator 
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for the hard probes, martini, to finally obtain a coupled 
simulation of both the soft and hard physics in heavy-ion 
collisions, creating an unprecedented theoretical tool for 
the study of the hot and dense phase of matter generated 
in heavy-ion collisions. 
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Appendix A: Kurganov-Tadmor Method 

As described in the main text, the Kurganov-Tadmor 
method is a MUSCL-type finite volume method in which 
the cell average of the density p around Xj is used instead 
of the value of the density at Xj. To do so, we first divide 
the space into equal intervals of the width Ax. If one 
integrates over the interval [xj — Ax/2,Xj + Ax/2], the 
conservation equation becomes 



dt 



•j + l/2 



dx p(x, t) = 

J(Xj-i/2,t) - J(x j+ i/ 2 ,t) 



(Al) 



where we introduced the notations £.,±1/2 = x j i Ax/2. 
Defining the cell average at Xj 0/S 



1 

Aa 



X 3 + l/2 



dx p(x, t) . 



(A2) 



'3,-1/2 

the above equation becomes 



d s , a _ J ( x j-i/2,t) - J(x j+1/2 ,t) 

Jt Pj W ~ Ax ' [A6) 

Using this exact equation, we can formally advance the 
time by At as 



p j (t + At) = p j (t) 



1 
~A~x 



t+At 



dt' (J(x j+1 / 2 ,f) - J{x J ^ 1 /2,t')) 



(A4) 



So far no approximation has been made. The main 
variables to calculate are the discrete average values 
Pj = f>j(t n ) for all Xj = xq + j Ax and at every time 
step t n = t + 11 At. Note that in Eq. (|A4|) . the cur- 
rents are evaluated at half-way points between Xj and 
its neighboring points Xj±\ . Therefore, even if Eq. (|A4[) 



is exact, it is not complete. One needs to know how to 
evaluate J{Xj±i/2,t') which in turn needs the value of 
p(xj±i/2,t'). Note further that this is not the cell av- 
erages but actual local values of p. Thus, the problem 
to solve now is how to approximate the current and the 
charge density at arbitrary x from the cell averages at 
discrete points Xj. 

A simple but effective solution is to approximate the 
local value with the average value and make a linear in- 
terpolation within each cell: 



p(x,t n )=J2[p] + (P*)? 



x 6(x j _ 1 /2 < x < x j+ i/ 2 ) (A5) 

where 0{xj-i/ 2 < x < Xj + i/ 2 ) is defined to be 1 when 
the condition is fulfilled and otherwise. This piecewise 
linear approximation is constructed in such a way that 
the amount of matter in the cell remains the same, that 
is, by construction 



dxp{x,t n ) = p" 



(A6) 



-1/2 



which, of course is a necessary condition for solving a 
conservation equation. The derivative (px)j 1 is a suitable 
approximation of d x p at Xj and t n constructed from the 
cell averages. It could be the backward slope 



Mi 



Ax 



or the forward slope 



Mi 



0+1 



(A7) 



(A8) 



or any combination of them. We will discuss the choice 
of derivatives in more detail shortly. For now we leave 
this choice open. 

There is, however, a potentially serious problem with 
this piecewise linear reconstruction. There are two ways 
to calculate the value of p(xj + i/ 2 ,t n ) = Pj +1 / 2 - One way 
is to calculate if from the left 



Pj+1/2 

or from the right 



left 



m + {p x )^Ax/2 { 



(A9) 



n' 6 



Ail 



J 3 + l 



( Px )"Ax/2. (A10) 



These two values in general do not coincide. Since we 
need these half-way values to calculate the currents at 
the cell boundaries, we need to find a way to deal with 
this discontinuity. 

A simple but effective solution to this problem was 
originally proposed by Nessyahu and Tadmor (74|: First, 
the discontinuity matters only for the current part be- 
cause of the spatial derivative. For the density part, 
the linear approximation is fine because only the time 
derivative is taken. For the density, we can integrate 
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over this discontinuity without any problem. Now, for a 
sufficiently small time interval At, the effect of the dis- 
continuity at xj + i/2 will not reach Xj and iCj+i. Hence, 
if one alternatively considers a cell defined by [x.j, Xj+i] 
instead of one defined by [^-1/21^+1/2]; then the cur- 
rents are calculated at Xj and Xj+i where p(x,t) is still 
smooth. Hence we can use p(x, t) to calculate both the 
charge density and the currents on the right hand side 
of Eq. (|A4|) provided that the average is now over the 
staggered cell [xj,Xj + i]. After this step, the discontinu- 
ities are located at Xj's instead of the half-way points. 
The next step in this approach is to repeat the same 
procedure now for [2^-1/2,2^+1/2] where this time the 
half-way points ^±1/2 are where the linear interpolation 
is smooth. This staggered grid approach is an effective 
method. However, the numerical viscosity term turns 
out to be 0((Aa;) 4 /At) which means that one still can- 
not take the At — > limit. 

Generalizing the Nessyahu-Tadmor method, Kurganov 
and Tadmor came up with a better solution to this prob- 
lem. Their main idea can be described as follows: The 
Nessyahu-Tadmor method relies on the smallness of At 
to guarantee that the influence of the discontinuities 
does not reach the mid-points of the (staggered) cells. 
This can be further improved if one uses one more piece 
of information, namely the maximum local propagation 
speed. The influence of the discontinuities at Xj±i/2 can 
travel no faster than the maximum propagation speed 
given by a = \dJ/dp\. Therefore, it makes sense to di- 
vide the space into two groups. One such group is given 
by the following set of intervals 

M"+i/2 = [x 3 +i/2-a] +1/2 At,x :j+1/2 +a] +1/2 At}, (All) 

where a" +i / 2 is the maximum propagation speed at 
Xj+1/2 an d time t n . The linear interpolation p(x,t) is 
possibly discontinuous in Mj l +i/2 as indicated by the 6- 
functions in Eq. (|A5|) . The other group is given by the 



'"two = 1 —ir- fd£p(£,tn) 1 — r - [ dt +/ Y 

J + l/2 on Af / ii srvs ' n > on Af I I 

M 3 + l/2^ 1 "'"■V1/2 Za 3 + l/2^ 1 Jin V 



where 2a" +1 y 2 At is the size of the interval. Using p(x, t) 
for the p integral and using the mid-point rule for the 
time integral, we obtain 

n" + p", , Ax - a/?, , /0 At 
<i/2=^ ±i + f^ 2 - ((*)? - (P^+i) 

-w-^(j(*?+i/2,+,t + &t/2) 

Z j+1/2 V 

-■7(*?+i/a,- ) * + At/2)) ) (A16) 

where we have defined 

x ]+i/2.± = X3+1/2 ± a] +1/2 At . (A17) 



set 

Xj = [xj-1/2 + a]_ 1/2 At, x j+1/2 - a] +1/2 At] , (A12) 

and in these intervals, p(x,t) is linear and smooth. The 
fact that we must have non-empty x™ gives us a condition 
on At 

Ax 

At < —n , (A13) 

a 3 + l/2 +a 3-l/2 

which is related to the CFL (Courant-Friedrichs-Lewy) 
condition. In fact, Ref. [tT| has a more severe CFL con- 
dition 

Ax 

At < , (A14) 

where a max is the maximum propagation speed in the 
whole region. 

The derivation of the Kurganov- Tadmor method now 
proceeds as follows. First, apply Eq. (|A4j) to M™ ±1 / 2 an d 
X™. Since we have divided the space into the //-set and 
the x-set, the currents are not being evaluated at the 
discontinuities at Xj±i/ 2 . Hence we can safely use p(x, t) 
to evaluate the right hand side of Eq. (|A4j) in both M"±i/2 
and Xj- I n this way, we get estimates of the density 
average in these intervals at the next time step. The 
intervals M™±i/2' s ancl Xj ' s are non-uniformly distributed. 
However, our starting point was the cell averaged values 
in the uniform grid. The next step is then to project 
the next-time values in this non-uniform grid of (J,j±i/ 2 's 
and Xj' s onto the original uniform grid [xj_i/ 2 ,Xj + i/ 2 ]. 
Finally, one takes the At — s- limit to get the semi- 
discrete equations. 

Application of Eq. (|A4|) to ^"±1/2 an d proceeds as 
follows. Within M™ +1 / 2 , the right hand side of Eq. (|A4[) 
is given by 



{xj+i/2 + a] +1/2 At, t) - J(x j+1/2 - a] +1/2 At, tfj (A15) 

I 

Similarly, for x™, 

w «+^p7-i( Px )7(a^ 1/2 -ay_ 1/2 )Ai 

(At/Ax) 

" l-(At/Ax)(aJ_ 1/2 + a« +1/2 ) 

x (j(x?+i/2,-,t + &t/2) 

-J(x";_ 1/2 , + ,t + At/2)Y (A18) 

At this point, w" +1 and w™+-y i 2 approximate the value 
of p at t n +i on a non-uniform grid. The next step is 
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to construct a piecewise linear function q(x,t n+ i) using Here 
w" +1 and w™+±i2 an d integrate over [£7-1/2,2^+1/2] to 
get the next cell average p™ +1 - For this purpose, the 
piecewise linear function g(x, is constructed by us- 

ing a linear approximation within p™ +1 , 2 and the con- 
stant approximation within Xj 

q(x,t n+1 ) = 

" {[wjfjj + (PxJTfVa^ - ^+1/2)] 



x 0(x e M?+ 1/a ) + g x")} (A19) 

The derivative appearing in the above approximation 
must also be calculated using w^^m an d w *j + ■ Again 
leaving what to use for the derivate for later discussions, 
integrating over the interval (xj-1/2, 2^+1/2) finally yields 
the cell average at the next time step 

n n+1 - lt, n+1 r, n —4- m n+1 n n — 
Pj ~ W j-l/2 a j-l/2 Aj; + W j+l/2 a j+l/2 Ax 



+ ( 1 -^K+i/2 + ^i/2))^; 

(a n l7 ,Ai) 2 
+ ( )™+i : ^ / _ 

_ / \n+l v J + l/2 ' 

lPxJj+1/2 2Ax 



n+l 

3 



(A20) 



Passing to the At — > limit, we get Kurganov and 
Tadmor's main result in the semi-discrete form 

IftW- ^V^ - (A21 » 



where 



if 



J±l/2 



•/(»j±l/2,+ )*) + ^0j±l/2 ,-,*) 



a j±l/2W 



(Pi± l /a,+ (*)-p 3 -±i/3 1 -(*)) (A22) 



£7+1/2,+ — Pj+i 



3+1 



Pj+1/2, 



Pj 



Ax 



(Pz)j 



(A23) 
(A24) 



and J{Xj±i/2 : ±) are evaluated with pj+1/2, ±- Any ex- 
plicit x dependence in pj+1/2, ± and J(xj±i/ 2 ,±) must be 
evaluated at Xj + i/ 2 - Note that all references to the in- 
termediate values have disappeared. 

One detail we need to take care of now is the choice 
of the spatial derivatives. Formally, the second-order ap- 
proximation 



(Ps) 



Pj+i ~ Pj-i 
2Aa; 



(A25) 



gives a better approximation than the first-order approx- 
imations Eq. (|A7[ IA8|) . But it is also known that when 
there is a stiff gradient, the second order expression tends 
to introduce spurious oscillations in the solution. To rem- 
edy this situation, one needs to use flux limitcrs which 
automatically switch the form of the numerical deriva- 
tive according to the stiffness of the slope. Kurganov 
and Tadmor chose the minmod limiter given by 



(p x )j = minmod 



where 



,Pi+i - Pj Pj+i - Pj-i Q P] ~ Pj-i 



2Aa 



Aa 



minmod (xi, x%, ■ 



muijjxj}, if Xj > Vj 
= i ma,Xj{xj}, if < Vj 
0, otherwise 



and 1 < 8 < 2 is a parameter that controls the amount 
of diffusion and the oscillatory behavior. This is also our 
choice with 6 = 1.1. 
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